The decay of turbulence in rotating flows 
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We present a parametric space study of the decay of turbulence in rotating flows combining direct 
numerical simulations, large eddy simulations, and phenomenological theory. Several cases are 
considered: (1) the effect of varying the characteristic scale of the initial conditions when compared 
with the size of the box, to mimic "bounded" and "unbounded" flows; (2) the effect of helicity 
(correlation between the velocity and vorticity); (3) the effect of Rossby and Reynolds numbers; 
and (4) the effect of anisotropy in the initial conditions. Initial conditions include the Taylor-Green 
vortex, the Arn'old-Beltrami-Childress flow, and random flows with large-scale energy spectrum 
proportional to k 4 . The decay laws obtained in the simulations for the energy, helicity, and enstrophy 
in each case can be explained with phenomenological arguments that consider separate decays 
for two-dimensional and three-dimensional modes, and that take into account the role of helicity 
and rotation in slowing down the energy decay. The time evolution of the energy spectrum and 
development of anisotropics in the simulations are also discussed. Finally, the effect of rotation and 
helicity in the skewness and kurtosis of the flow is considered. 



I. INTRODUCTION 

Nature presents several examples of rotating flows. Ro- 
tation influences large-scale motions in the Earth's atmo- 
sphere and oceans, as well as convective regions of the sun 
and stars. Rotation is also important in many industrial 
flows, such as turbo machinery, rotor-craft, and rotating 
channels. In a rotating system, the Coriolis force, lin- 
ear in the velocity, modifies the flow nonlinear dynamics 
when strong enough. In its presence, the Navier-Stokcs 
equation becomes a multi-scale problem with a "slow" 
time scale tl ~ L/U associated with the eddies at a 
characteristic scale L (U is a characteristic velocity), and 
a "fast" time scale tq ~ I/O ~ t^Ro associated with 
inertial waves. The dimensionless Rossby number Ro is 
the ratio of advection to Coriolis forces, and measures 
the influence of rotation upon the nonlinear dynamics of 
the system (decreasing as rotation becomes dominant). 

Resonant wave theory provides a framework to study 
the effect of rapid rotation in turbulence 0-01 ■ The sep- 
aration between the fast and slow time scales results in a 
selection of the resonant triadic interactions as the ones 
responsible for the energy transfer among scales. As a 
result, energy transfer and dissipation is substantially 
decreased in the presence of strong rotation Q. The 
resonant condition is also responsible for the transfer of 
energy towards two-dimensional (2D) slow modes, driv- 
ing the flow to a quasi-2D state [3, 0] (this result is of- 
ten referred to in the literature as the "dynamic Taylor- 
Proudman theorem", see e.g., @). The development of 
anisotropy and reduction of the energy transfer and dis- 
sipation rates has been verified in numerical simulations 
[JClI and experiments (HEl. 

Similar arguments (see, e.g., [HI) indicate that in 
the limit of fast rotation (small Rossby number) the 
slow 2D modes decouple from the remaining fast three- 



dimensional (3D) modes, and evolve under their own au- 
tonomous dynamics. Moreover, in that limit the av- 
eraged equation for the slow modes splits into a 2D 
Navier-Stokes equation for the vertically-averaged hor- 
izontal velocity and a passive scalar equation for the 
vertically-averaged vertical velocity. Although simula- 
tions of forced rotating flows @ and of ideal truncated 
rotating flows fl5l.[l6| using periodic boundary conditions 
show good agreement with these predictions for small 
enough Rossby numbers, for long times the decoupling 
of slow and fast modes seems to break down. Also, some 
authors [TtJ argue that in unbounded domains no decou- 
pling is achievable even for Ro = 0. 

Of particular importance for many geophysical and as- 
trophysical problems is how turbulence decays in time. 
The problem is also important for laboratory experi- 
ments, as it provides, e.g., one way to measure changes in 
the energy dissipation rate associated with the presence 
of rotation. Even in the absence of rotation, the decay of 
isotropic turbulence proves difficult to tackle because of 
the different decay laws obtained depending on boundary 
and initial conditions. As an example, for bounded flows 
(i.e., flows for which the initial characteristic length is 
close to the size of the vessel) the energy decays as ~ t~ 2 
(see, e.g., [U, EH EH)- For unbounded flows (i.e., flows 
in an infinite domain, or in practice, flows for which the 
initial characteristic length is much smaller than the size 
of the vessel) a - t~ w / 7 @, UK or ~ £" 6 / 5 [H HI H 
decay law is observed depending on whether the initial 
energy spectrum at large scales behaves as ~ k A or ~ k 2 , 
respectively. 

In the presence of rotation, the decay of turbulence 
becomes substantially richer, with decay rates depend- 
ing not only on whether turbulence is bounded or un- 
bounded and on the initial spectrum at large scales, but 
also, e.g., on the strength of background rotation (see 
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|24|). A detailed experimental study of this dependence 
can be found in Refs. [H, [H|, where the authors studied 
the energy decay of grid-generated turbulence in a ro- 
tating tank using particle image vclocimetry, and found 
different decay laws depending upon the rate of rotation 
and the saturation (or not) of the characteristic size of the 
largest eddies. For large Rossby number they reported 
a decay ~ t a with exponent a w —1.1 for non rotating 
turbulence (consistent with the value of —6/5 predicted 
in [22|, [23[ for ~ k 2 unbounded turbulence), which later 
turned to a w — 2 after the largest eddies grew to the 
experiment size. For small Rossby number this energy 
decay rate became smaller saturating at « —1. Simi- 
lar results were reported in simulations in [Toj . were a 
decrease from rss —10/7 to —5/7 was observed as rota- 
tion was increased. These results are consistent with the 
reduction of the energy transfer discussed above. They 
also reported a steeper energy spectrum together with 
positive vorticity skewness for the rotating flows, and 
anisotropic growth of integral scales (see also 0, [25| ) . 

The decay rate of rotating turbulence also seems to 
depend on the helicity content of the flow. Helicity (the 
alignment between velocity and vorticity) is an ideal in- 
variant of the equations of motion (a quantity conserved 
in the inviscid limit) with intriguing properties. In the 
absence of rotation, the presence of helicity does not 
modify the energy decay rate [ljl [2|| HtJ nor the dissipa- 
tion rate (2^. However, in the presence of rotation [22} 
reported a further decrease of the energy transfer when 
both rotation and helicity were present, and [l9[ showed 
that the decay rate of bounded rotating flows changes 
drastically depending on whether helicity is present or 
not. 

In this paper we conduct a detailed study of param- 
eter space of rotating helical flows, taking into account 
(1) the effect of varying the characteristic scale of the 
initial conditions when compared with the size of the 
box, (2) the effect of helicity, (3) the effect of Rossby 
and Reynolds numbers, and (4) the effect of anisotropy 
in the initial conditions. The numerical study uses a 
two pronged approach combining direct numerical simu- 
lations (DNS) and large eddy simulations (LES). Several 
initial conditions are considered, although when the char- 
acteristic initial scale of the flow is smaller than the size 
of the domain, we focus only on the case with large scale 
initial energy spectrum ~ fc 4 . The different decay laws 
obtained (which in some cases coincide with previous ex- 
perimental or numerical observations, while in others are 
new) are explained using phcnomenological arguments, 
and we classify the results depending on the relevant ef- 
fects on each case. 

After studying the decay laws, we study the evolution 
of anisotropy, of skewness and kurtosis, and the forma- 
tion of columnar structures in the flow. We consider 
how helicity affects the evolution of skewness and kur- 
tosis, and associate peaks observed in the time evolution 
of these quantities with the dynamics of the columnar 
structures. 



The structure of the paper is as follows. In Sec. IPTl we 
introduce the equations, describe the DNS and LES, and 
give details of the initial conditions and different param- 
eters used. Section [TTT] presents phcnomenological argu- 
ments to obtain decay laws in turbulent flows with and 
without rotation. The phenomenological predictions are 
then compared with the numerical results in Sec. IIV1 
which presents the energy, helicity, and enstrophy decay 
in all runs. The spectral evolution and development of 
anisotropy in the flows is discussed at the end of that 
section. The effect of initial anisotropy in the decay is 
considered in Sec.|V] A statistical analysis including evo- 
lution of skewness and kurtosis is presented in Sec. I VII 
Finally, Sec. IVIII gives our conclusions. 



II. EQUATIONS AND MODELS 

A. Equations 

The evolution of an incompressible fluid in a rotating 
frame is described by the Navier-Stokes equation with 
the Coriolis force, 

9fU + wxu + 2nxu = -VP+i/V 2 u, (1) 

and the incompressibility condition, 

V-u = 0, (2) 

where u is the velocity field, u> = V x u is the vorticity, 
the centrifugal term is included in the total pressure per 
unit of mass V, and v is the kinematic viscosity (uni- 
form density is assumed). The rotation axis is in the z 
direction, so fi = Oz, where Q is the rotation frequency. 

As mentioned in the introduction, these equations are 
solved numerically using two different methods: DNS, 
and LES using a dynamical subgrid-scale spectral model 
of rotating turbulence that also takes into account the 
helicity cascade. All simulations were performed in a 
three dimensional periodic box of length 27r, using differ- 
ent spatial resolutions ranging from 96 3 grid points for 
the lowest resolution LES runs up to 512 3 for the highest 
resolution DNS. 



B. Models 

In a DNS, all spatial and time scales (up to the dis- 
sipation scale) are explicitly resolved. The simulations 
were performed using a parallelized pseudo-spectral code 
[29l [30| with the two-third rule for dcaliasing. As a re- 
sult, the maximum wave number resolved in the DNS is 
kmax — N/3 where N is the linear resolution; to properly 
resolve the dissipative scales the condition k^/kmax < 1 
must be satisfied during all simulations, where k v is the 
dissipation wave number. In practice, this condition is 
more stringent if reliable data about velocity gradients 
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TABLE I: Parameters used in the simulations: kinematic viscosity v, rotation rate Q, Reynolds number Re, Rossby number 
Ro, micro- Rossby number Ro" , initial relative helicity h, relative helicity at the time of maximum dissipation h* , and time 
of maximum dissipation t* . The values of Re, Ro, and Ro" are always given at t* . The last column succinctly describes the 
initial energy spectrum E(k): the power law followed by the spectrum, the range of scales where this power law is satisfied, 
and the flow (TG for Taylor-Green, ABC for Arn'old-Beltrami-Childress, and RND for random). 



Run 


V 


Q. 


Re 


Ro 


Ro" 


h 


h* 


t* 


Initial E(k) 


D256-1 


1.5 x l(T a 





450 


CO 


CO 





9 x 10" 1U 


1.26 


fc" 4 (4-14) TG 
k~ 4 (4-14) TG 


D256-2 


1.5 x 10" 3 


4 


550 


0.12 


1.28 





-1 x 10" 8 


1.06 


D256H-1 


1.5 x 1(T 3 





600 


oo 


CO 


0.95 


0.34 


2.28 


AT 4 (4-14) ABC 
k~ 4 (4-14) ABC 
k~ 4 (4-14) TG 


D256H-2 


1.5 x 10" 3 


4 


830 


0.08 


0.80 


0.95 


0.65 


2.25 


D512-1 


7 x 10~ 4 


4 


1100 


0.12 


1.82 





7 x 10~ 9 


0.88 


D512-2 


8.5 x 1CT 4 





420 


CO 


CO 


8 x 10~ 5 


8 x 10" 4 


0.60 


fe 4 (1-14) RND 
k 4 (1-14) RND 
k~ 4 (4-14) ABC 


D512-3 


8.5 x itr 4 


10 


450 


0.10 


0.95 


4 x 10~ 3 


4 x 10~ 3 


0.70 


D512H-1 


7 x 10~ 4 


4 


1750 


0.08 


1.15 


0.95 


0.44 


1.70 


D512H-2 


8 x 1(T 4 





440 


CO 


CO 


0.90 


0.38 


0.94 


k 4 (1-14) RND 


D512H-3 


8 x 10" 4 


10 


530 


0.07 


0.70 


0.99 


0.5 


1.50 


k 4 (1-14) ABC 




o.O X 1U 


u 


^ pin 


OO 


X. 


n 


n no 


n 


K ^1-14^ it IN U 


L96-2 


8.5 x itr 4 


2 


540 


0.42 


2.90 


-0.03 


-0.02 


0.30 


fe 4 (1-14) RND 


L96-3 


8.5 x itr 4 


4 


540 


0.21 


1.45 


-0.03 


-0.02 


0.30 


k 4 (1-14) RND 


L96-4 


8.5 x itr 4 


6 


550 


0.14 


0.95 


-0.03 


-0.02 


0.30 


k 4 (1-14) RND 


L96-5 


8.5 x itr 4 


8 


550 


0.11 


0.73 


-0.03 


-0.02 


0.30 


k 4 (1-14) RND 


L96-6 


8.5 x itr 4 


10 


530 


0.08 


0.65 


0.03 


0.02 


0.30 


k 4 (1-14) RND 


L96H-1 


8 x 1(T 4 





500 


CO 


CO 


0.90 


0.51 


0.70 


k 4 (1-14) RND 


L96H-2 


8.5 x itr 4 


10 


540 


0.08 


0.63 


0.90 


0.70 


0.70 


k 4 (1-14) RND 


L96H-3 


8.5 x 10~ 4 


10 


490 


0.08 


0.60 


0.99 


0.80 


1.15 


k 4 (1-14) ABC 


L192-1 


2 x 10~ 4 





1200 


CO 


CO 


-7 x 10" 3 


-6 x 10" 3 


0.10 


k 4 (1-30) RND 


L192-2 


2 x 1(T 4 


10 


1100 


0.22 


1.65 


-7 x 10~ 3 


-6 x 10~ 3 


0.13 


k 4 (1-30) RND 


L192H-1 


2 x 1(T 4 





950 


CO 


CO 


0.90 


0.60 


0.30 


k 4 (1-30) RND 


L192H-2 


2 x 1(T 4 


10 


1000 


0.20 


1.60 


0.94 


0.71 


0.38 


k 4 (1-30) ABC 


L192HA-1 


2 x 1(T 4 


10 


1200 


0.16 


1.40 


0.90 


0.56 


0.50 


k 4 (1-25) RND 


L192HA-2 


2 x 1(T 4 


10 


1300 


0.14 


1.35 


0.90 


0.59 


0.46 


k 4 (1-25) RND 


L192HA-3 


2 x 10~ 4 


10 


1300 


0.15 


1.35 


0.90 


0.58 


0.45 


k 4 (1-25) RND 



and high-order statistics of the flow are needed (see, e.g., 
31 1, where they indicate k^/kmax < 0.5 for such studies). 

The dissipation wave number as a function of time was 
computed for all simulations in two different ways: as the 
Kolmogorov dissipation wave number for isotropic and 
homogeneous turbulence k^ — (e/i^ 3 ) 1 / 4 = ((o; 2 )/^ 2 ) 1 / 4 
(where e is the energy dissipation rate and (w 2 ) the mean 
square vorticity), and as the wave number where the en- 
strophy spectrum peaks. The Kolmogorov dissipation 
wave number was found to be always larger than the 
wave number where dissipation peaks, and in the follow- 
ing we therefore only consider the Kolmogorov scale as 
it gives a more stringent condition on the resolution. 

In all DNS discussed below, the ratio k n jk max was 
< 0.7 at the time of maximum dissipation (t k 1 to 
tw3 depending on the simulation), 0.2 to 0.5 at t ~ 10 
(when the self-similar decay starts in most runs), and 
monotonously decreases to values between 0.05 to 0.2 at 
t ss 100. The spatial resolutions used were 256 3 and 512 3 
grid points, and an explicit second-order Runge-Kutta 
method was used to evolve in time, with a Courant- 
Friedrichs-Levy (CFL) number smaller than one. 

In the LES approach, only the large scales are explic- 
itly resolved, while the statistical impact on the resolved 



scales of scales smaller than a cut-off scale is modeled 
with simplified equations. To this end we use the spectral 
model derived in (32| for isotropic helical and non-helical 
turbulence, and its extension to the rotating case [33| . 
The model is based on the eddy damped quasi-normal 
Markovian (EDQNM) closure to compute eddy viscos- 
ity and eddy noise, and assumes unresolved scales (scales 
smaller than the cut-off) are isotropic. Both eddy viscos- 
ity and noise are computed considering the contribution 
from the energy and the helicity spectra (see, e.g., [34| 
for another subgrid model that takes into account the 
effect of helicity). The model adapts dynamically to the 
inertial indices of the resolved energy and helicity spec- 
tra, and as a result it is well suited to study rotating 
turbulence for which the scaling laws are not well known 
and may depend on the Rossby number. For a validation 
of the LES against DNS results, the reader is referred to 

The subgrid model starts by applying a spectral fil- 
ter to the equations; this operation consists in truncat- 
ing all velocity components at wave vectors k such that 
|k| = k > k c , where k c is the cut-off wave number. One 
then models the transfer between the large (resolved) 
scales and the small (subgrid unresolved) scales of the 
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flow by adding eddy viscosity and eddy noise to the equa- 
tions for the resolved scales. These are obtained solving 
the EDQNM equations for estimated energy and helicity 
spectra in the subgrid range. To this end, an interme- 
diate range is defined, lying between k' c and k c (in most 
cases k' c = fc c /3), where the energy spectrum is assumed 
to present a power-law behavior possibly followed by an 
exponential decrease. As an example, for the energy the 
following expression is used: 



E(k,t) = E k~ 



k r . ^ k kf< 



(3) 



The coefficients ag, Se and Eq are computed at each 
time step doing a mean square fit of the resolved energy 
spectrum. The spectrum is extrapolated to the unre- 
solved scales using these coefficients, and the EDQNM 
equations are solved. Then one solves the Navier-Stokcs 
equation (JTJ with an extra term on the r.h.s. which in 
spectral space takes the form 



-v{k\k c ,t) fc 2 u(k,t), 



(4) 



where u(k, t) is the Fourier transform of the velocity 
field u(x, t), —k 2 is the Laplacian in Fourier space, and 
v (k\k c , t) is an eddy viscosity proportional to the ratio of 
the so-called absorption terms in EDQNM to the energy 
(and helicity) spectrum. Eddy noise is added in a similar 
manner (for more details, see [32l l33j). 

LES runs using this model have a resolution of cither 
96 3 or 192 3 grid points. A pseudo-spectral method is also 
used, but without dcaliasing, resulting in the maximum 
wave number k max = N/2. As in the DNS, an explicit 
second-order Rungc-Kutta method is used to evolve in 
time. 

Parameters for all sets of runs are listed in Table Q] 
DNS runs are labeled with a D, followed by the linear 
resolution, a letter "H" if the run has helicity, a letter "A" 
if the initial energy spectrum is anisotropic, and the run 
number. LES runs start with an L, followed by numbers 
and letters using the same convention as in the DNS. 



C. Initial conditions and definitions 

As mentioned in the introduction, we are interested 
in the decay laws obtained in the system depending on 
properties of the initial conditions and the amount of 
rotation. In particular, we will vary the initial amount of 
helicity, the initial energy containing scale (with respect 
to the largest available scale in the box), the shape of 
the energy spectrum, and the strength of turbulence and 
of rotation as controlled by the Reynolds and Rossby 
numbers respectively. 

Helicity is an ideal invariant of the Navier-Stokes equa- 
tion which measures the alignment between velocity and 
vorticity. If zero, the initial conditions are mirror- 
symmetric, so it also measures the departure from a 
mirror-symmetric state. We define the net helicity as 



where the brackets denote spatial average. We also define 
the relative helicity as 



h 



H 



(0) 



MM)' 

which is bounded between —1 and 1 and can be inter- 
preted as the mean cosine of the angle between the ve- 
locity and the vorticity. 

To control the net amount of relative helicity in the 
initial conditions we consider three different flows: a su- 
perposition of Taylor-Green (TG) vortices (HJ, 

U sui(kox) cos(fcoy) cos(koz)x — 
U cos{k x) sin(fc y) cos(koz)y, (7) 



Utg 



a superposition of Arn'old-Beltrami-Childrcss (ABC) 
flows 13611. 



vlabc = [Bcos(k y) + Asm(k z)]x-\ 
[C cos(k z) + A sin(fc a;)] y -\ 
[A cos(fc x) + B sin(fc 2/)] z, 



(8) 



(with A = 0.9, B = 1.1, and C = 1), and a superposition 
of Fourier modes with random phases (RND) in which we 
use the algorithm described in [37j to control the relative 
helicity. In each case, the flows are superposed in a range 
of wave numbers as described in Table U and with global 
amplitudes for each wave number to give the desired slope 
in the initial energy spectrum. 

The TG vortex is non-helical (h = 0) and has no en- 
ergy in the k z = modes, whose amplification in the 
rotating cases (see below) can thus be attributed only to 
a bidimensionalization process. The TG vortex was orig- 
inally motivated as an initial condition which leads to 
rapid development of small spatial scales, and also mim- 
ics the von Karman flow between two counter-rotating 
disks used in several experiments [38[. The ABC flow is 
an eigenfunction of the curl operator and as a result has 
maximum helicity (h = ±1 depending on the sign of kg, 
when only one value of ko is excited), whereas the RND 
flow allows us to tune the amount of initial relative he- 
licity between —1 and 1 as well as the initial anisotropy. 

When generating the flows, two different initial energy 
spectra were considered. To study initial conditions with 
characteristic length close to the size of the computa- 
tional domain, a spectrum E(k) ~ fc~ 4 for k G [4, 14] 
(followed by exponential decay) was imposed. To study 
initial conditions with length smaller than the domain 
size, we imposed a spectrum E(k) ~ k 4 for k € [1,14], 
[1,25], or [1,30] (also followed by exponential decay). In 
the latter case, the characteristic length can grow in time 
as the spectrum peaks around k = 14, 25, or 30. This 
allows us to mimic (at least for a finite time before the 
characteristic length reaches the domain size) the decay 
of unbounded flows. The characteristic length will be 
associated in the following with the flow integral scale, 
which is defined as 



H = (u-u), 



(5) 
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E k E(k) ' 



(9) 
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where E(k) is the isotropic energy spectrum. 

Simulations in Table U are also characterized by differ- 
ent Reynolds and Rossby numbers. The Reynolds and 
Rossby numbers in the Table are defined as 



and 



Re 



Ro 



UL 
v 

U 



(ii) 



respectively Of importance is also the micro-scale 
Rossby number (see e.g., 0) 



Rn"> - w 

R ° -2Q' 



(12) 



which can be interpreted as the ratio of the convective 
to the Coriolis acceleration at the Taylor scale. The 
Rossby number Ro must be small enough for rotation 
to affect the turbulence, while the micro-Rossby number 
Ro u must be larger than one for scrambling effects of in- 
ertial waves not to completely damp the nonlinear term, 
which would lead to pure exponential viscous energy de- 
cay Q . In all runs in Table |fl Ro and i?o w are one order 
of magnitude apart at the time of maximum enstrophy 
t* , and this interval is roughly preserved throughout the 
simulations. 

Here and in the following, the isotropic energy spec- 
trum is defined by averaging in Fourier space over spher- 
ical shells 



E{k,t) 



u*(k,t).u(k,t), (13) 



fe<|k|<fe+l 



where u(k, t) is the Fourier transform of the velocity field, 
and the asterisk denotes complex conjugate. Other two 
spectra can also be used to characterize anisotropy. 

On the one hand, the so-called "reduced" energy spec- 
tra E(kj_) and E(k\\) arc defined averaging in Fourier 
space over cylinders and planes respectively. More specif- 
ically, the reduced energy spectra as a function of wave 
numbers k± with = (k x ,k y ,0), and kit with k|| = 
(0,0, k z ), are defined by computing the sum above over 
all modes in the cylindrical shells k± < |kj_| < k± + 1 
and over planes fcii < |kj| | < fcy + 1 respectively (isotropic 
and reduced spectra for the helicity are defined in the 
same way). From the reduced spectra, perpendicular and 
parallel integral scales can be defined; e.g., for the per- 
pendicular direction, 
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On the other hand, more information of the spectral 
anisotropy can be obtained studying the axisymmetric 
energy spectrum e(fc||,fcj_) (see, e.g., [H, HI). Assuming 
the flow is axisymmetric, the three-dimensional spectrum 
can be integrated around the axis of rotation to obtain a 



spectrum that depends only on k» and k±, which relates 
to the reduced energy spectra as follows: 



(10) and 



s (*l|)=£ e (*l|.*-i). 
E(k ± ) = J2 e ( k h k ^)- 



(15) 



(16) 



III. TIME EVOLUTION - PHENOMENOLOGY 

We present now phenomenological arguments that will 
become handy to understand the different decay rates 
that are observed in our simulations as well as in previous 
studies. Some of the arguments are well known, while 
others are new, and we quote previous derivations when 
needed. 



A. Non- rotating flows 

1. Bounded 

From the energy balance equation 
dE 

~dt ~ 6 



(17) 



where e is the energy dissipation rate, Kolmogorov phe- 
nomenology leads to 



dE 
~dt 



£3/2 

— T 



(18) 



where E = E(t) ~ kE(k) and L is an energy containing 
length scale. For bounded flows where L ~ Lq (Lq is 
the size of the box), Eq. ([Tg]) becomes dE/dt ~ E 3 / 2 /L Q , 
resulting in the self-similar decay 0, HI, Eflj 



E(t) - t' 



(19) 



2. Unbounded 



In unbounded flows, a similarity solution of Eq. (|18[) 
requires some knowledge of the behavior of the energy 
containing scale L, which is in turn related to the evo- 
lution of E(k) for low wave numbers. In the case of an 
initial large scale spectrum ~ k\ the quasi-invariance of 
Loitsyanski's integral / (see [U, |4l[) leads, on dimen- 
sional grounds, to / ~ L 5 E, and replacing in Eq. (|18l) we 



(14) get Kolmogorov's result [20| 



E(t) - t 



-10/7 



(20) 



A different decay law is obtained if an initial ~ k 2 spec- 
trum is assumed for low wave numbers (22l . [23j . In the 
following, we will consider only the bounded or the ~ fc 4 
unbounded cases. 
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FIG. 1: (a) Energy decay for non-rotating unbounded runs. 
Non-helical runs D512-2 (solid), L96-1 (dashed), L192-1 
(dash-dotted), and helical runs D512H-2 (solid, thick), L96H- 
1 (dashed, thick), and L192H-1 (dash-dotted, thick) are 
shown. A —10/7 slope is shown as a reference. The inset 
shows the energy decay for non-rotating bounded runs D256- 
1 (solid), and D256H-1 (solid, thick), (b) Enstrophy decay 
for the same runs, with a —17/7 slope shown as a reference. 
The inset shows the enstrophy decay in the bounded runs, (c) 
Helicity decay in the unbounded helical runs of (a); the inset 
shows the helicity decay for the bounded helical run D256H-1. 



FIG. 2: Energy decay for different values of Q from to 10 
for unbounded, non-helical runs L96-1, L96-2, L96-3, L96-4, 
L96-5, and L96-6. The decay becomes slower with increasing 
rotation rate. We also show t~ 10 ^ 7 and t~ 5 ^ 7 laws as refer- 
ences. 



B. Rotating flows: isotropic arguments 

1. Bounded 

In the case of solid-body rotation without net helicity, 
a spectra E(k) ~ e 1 / 2 ^! 1 / 2 /^ -2 is often assumed at small 
scales (i-e., wave numbers larger than the integral wave 
number) 0, I421 - I451 ] . Replacing in the balance equation, 
this spectrum leads to 



dE 
~dt 



E 2 



(21) 



For bounded flows L ~ L a and we get [ij, Ej, |46| 

E(t)~i _1 . (22) 

In helical rotating flows the small-scale energy spec- 
trum takes a different form. The direct transfer is dom- 
inated by the helicity cascade. In this case we can write 
the helicity flux as 6 ~ hi/{Q,rf) where hi is the helicity 
at the scale I, and t\ the eddy turnover time [ll|. Con- 
stancy of S leads to small scale spectra E(k) ~ k~ n and 
H(k) ~ fc™ -4 , where n = 5/2 obtains for the case of max- 
imum helicity [ll|. Further use of dimensional analysis 
leads to E(k) - e 1 / 4 ^ 4 ^ 5 ? 2 for the energy spectrum 
in terms of the energy dissipation rate, and replacing in 
the balance equation we get 



dE 
~dt 



E 4 



(23) 



For L ~ L then hj 



E{t) - t 



'1/3 



(24) 
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2. Unbounded 



For non-helical unbounded flows with E{k) ~ k at 
small wave numbers we can again make use of the con- 
stancy of I in Eq. (|2"Tj) . leading to [25| 



-5/7 



(25) 



For helical flows, assuming / remains constant in 
Eq. we obtain [47 1 



£(t) - f 



-5/21 



(26) 
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FIG. 3: (a) Energy decay for rotating unbounded runs (f2 = 
10). Non-helical runs D512-3 (solid), L96-6 (dashed), L192-2 
(dash-dotted), and helical runs D512H-3 (solid, thick), L96H- 
3 (dashed, thick), L96H-2 (dash-tripe-dotted, thick), and 
L192H-2 (dash-dotted, thick) are shown. At late times, the 
non-helical runs decay slightly faster than t _5//7 , while the he- 
lical runs are close to a —1/3 decay. The inset shows bounded 
non-helical runs D256-2 (solid) and D512-1 (dashed), and he- 
lical runs D256H-2 (solid, thick) and D512H-1 (dashed, thick), 
(b) Enstrophy decay for the same runs, with a —12/7 slope 
shown as a reference. The inset shows the enstrophy decay 
in the bounded runs, (c) Helicity decay in the unbounded 
helical runs; bounded helical runs are shown in the inset. 



C. Rotating flows: anisotropic arguments 

The decay laws obtained for rotating flows in Eqs. (|22[) . 
(|24p . and (I25[). have been reported in experiments or in 
simulations jl2l Il3l. [l9l |2"H . However, the analysis above 
is based on the isotropic energy spectrum and on the 
quasi-invariance of the isotropic Loitsyanski integral. For 
an anisotropic flow, other quantities are expected to be 
quasi-invariants during the decay instead [4l|, HH, [49[ . 

Rotating flows tend to become quasi-2D, and the as- 
sumption of an axisymmetric energy spectrum seems nat- 
ural considering the symmetries of the problem. If there 
is no dependence on wave numbers on the parallel direc- 
tion, the energy spectrum for small values of k± can be 
expanded as (see, e.g., fofjj ) 

E{k±) ^ Lk^ 1 +Kk ± +I 2D k 3 ± + --- (27) 

We will be interested in the following coefficients: 



and 




where (u • u'} is the two-point correlation function for 
spatial increments r perpendicular to the rotation axis. 
If the correlation function decays fast enough for large 
values of r [H|, these quantities can be expected to be 
quasi-invariants during the decay, respectively for initial 
large-scale energy spectra ~ k± and ~ kj_, in the same 
way I is quasi-conserved during the decay of isotropic 
flows with an initial large-scale ~ k A energy spectrum. A 
detailed proof of the conservation of K for rotating flows 
can be found in ■ 49;] ; it is a direct consequence of the con- 
servation of linear momentum in the direction parallel to 
the rotation axis. It is worth pointing out that these 
quantities were also shown to be conserved in other sys- 
tems: proofs of the conservation of K and I2D for quasi- 
geostrophic flows can be found in [5(| . In practice, these 
quantities are only approximately conserved in numerical 
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FIG. 4: Energy decay for Ezd (solid) and E 2 o (dashed) for 
runs with rotation: (a) Unbounded non-helical D512-3 (thin) 
and L96-6 (thick); E 3 d ~ t~ 10/7 and E 2D ~ t~ 2/3 decays are 
indicated, (b) Unbounded helical with random initial con- 
ditions L96H-2; E 2 d is close to t -1,/2 . (c) Unbounded heli- 
cal with ABC initial conditions D512H-3 (thin) and L96H-3 
(thick); E 2 d is close to i~ 1/3 . (d) Bounded helical with ABC 
initial conditions D512H-1. 



time 



FIG. 5: Evolution of 7/7(0) (thick lines) and of I 2D /I 2 d{0) 
(thin lines) for runs D512-3 (solid) and L96-6 (dashed). While 
in both runs I 2 d mantains an approximately constant value, 
I growths monotonically and during the self-similar energy 
decay increases by one order of magnitude. 



simulations, see e.g., the approximate constancy of I 2 d 
and K reported for rotating flows in j47j . 

As per virtue of the decay the Rossby number decreases 
with time, we will further assume for our phcnomcnolog- 
ical analysis that 2D and 3D modes are only weakly cou- 
pled, and write equations for the energy in the 2D modes, 
E-2D ■ In the non- helical case, if K remains approximately 
constant with K ~ E 2 dE\L^ (where LqII is the size of 
the box in the direction parallel to Q), then Eq. (j2l"j) for 
the 2D modes becomes 



dE 2D E$ D L 



'II 



dt KQ 
which leads to a decay 

E 2D (t)~r 



1/2 



(30) 



(31) 



Alternatively, constancy of I2D ~ E 2 dL <i_L \\ in Eq. ^1 
for the 2D modes leads to 



dE' 



-ID 



dt 



p 5/2 r l/2 
^2D ^0\\ 

1 2D " 



and 



E 2D {t)~t 



-2/3 



(32) 



(33) 



The same arguments can be extended to the helical 
rotating case using Eq. (|23|) . If constancy of K is assumed 
we get 



dE- 



2D 



dt 



E 2D L 0\\ 



and 



E(t) - t 



-1/6 



(34) 



(35) 
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Finally, constancy of I2D leads to 

p,ll/2 r 3/2 

and 

E(t)~t- 2 / 9 . (37) 

These decay laws will be important to analyze the evo- 
lution of the energy in the simulations discussed in the 
next section. 



D. Enstrophy decay 

From any of the previous energy decay laws, one can 
also compute laws for the enstrophy decay = (w 2 ) /2 
using the isotropic energy balance equation and replacing 
e = vQ.(f), which results in f2(f) = v~ 1 dE/dt. From this 
equation, for every solution for which the energy decays 
as E(t) ~ t a , the enstrophy decay results 

n(t)~t a_1 . (38) 

Although rotating flows are anisotropic, the enstrophy 
is predominantly a small-scale magnitude and we will see 
that this isotropic argument gives good agreement with 
the numerical results for rotating and non-rotating flows. 
Since helicity is related with the energy and the enstro- 
phy only through a Schwartz inequality, no simple decay 
laws can be derived in its case using these phenomeno- 
logical arguments. 

IV. TIME EVOLUTION - NUMERICAL 
RESULTS 

We present here the results for the energy, enstrophy, 
and helicity decay obtained in the numerical simulations 
listed in Table HI classifying them as rotating or non- 
rotating, bounded or unbounded (in the sense that the 
initial integral scale is smaller than the size of the box), 
and helical or non-helical. 

Concerning the terminology of "bounded" and "un- 
bounded" used to describe the numerical simulations, it 
is important to note that confinement effects in a ro- 
tating flow go beyond a saturation of the integral scale 
when it grows to the box size. Confinement also selects 
a discrete set of inertial waves which are normal modes 
of the domain, and boundaries can introduce dissipation 
through Ekman layers. The latter effect is not present in 
our numerical simulations with periodic boundary con- 
ditions. Finally, it was shown in [5l[ (see also [52|) that 
the small number of Fourier modes available in the shells 
with wave number k « 1 gives rise to poor represen- 
tation of isotropy and of the integral scale in runs for 
which the integral scale approaches 1/5 of the box size. 



As a result, the "bounded" runs are here only briefly con- 
sidered to study the time evolution of global quantities 
(energy, enstrophy, and helicity), and to compare with 
the prediction obtained in the corresponding cases in the 
phcnomcnological analysis. 

A. Non-rotating flows 

Numerical results for non-rotating, bounded and un- 
bounded flows are shown in Fig. [1] In the unbounded 
case (runs with an initial energy spectrum ~ fc 4 peak- 
ing at k = 14 in the DNS and 96 3 LES, and peaking 
at k = 30 in the 192 3 LES), the runs show a decay for 
the energy close to ~ t -10 / 7 independently of the pres- 
ence of helicity or not (note the runs also span a range of 
Reynolds numbers from Re ps 420 to 1200). The decay 
is consistent with the prediction given by Eq. (|2T)]) for an 
initial ~ fc 4 energy spectrum. 

The enstrophy decay is also consistent with this law, 
as expressed by Eq. (f3"B")) . decaying close to ~ t~ 17 / 7 in 
all cases. In the absence of rotation, helicity only delays 
the onset of the self-similar decay by retarding the time 
when the maximum of enstrophy takes place, as already 
reported in (2?| and [53|. This is more clearly seen in 
the DNS; see, e.g., the time of the peak of enstrophy for 
runs D512-3 and D512H-3 in Fig. QTb). Finally, in the 
helical runs, helicity seems to decay as the enstrophy, just 
slightly slower than the ~ t~ 17 ^ 7 law. 

Similar results are observed for bounded flows, i.e., for 
initial conditions with the initial energy containing scale 
close to the size of the box (runs with a ~ k~ A spectrum 
from k — 4 to 14, peaking at k = 4). In this case, all 
runs are consistent with a ~ t~ 2 decay for the energy 
(see the insets of Fig. QJ in agreement with Eq. (|T9|) . and 
a decay for the enstrophy close to ~ t~ 3 in agreement 
with Eq. (f3"5| . In the helical runs, helicity decays again 
slightly slower than the enstrophy, but close to the ~ t~ 3 
power law. 

B. Rotating flows 

1. Global quantities 

As rotation is increased, the simulations show a shal- 
lower power law in the energy decay. As an illustration, 
Fig. [5] shows the energy decay rate in simulations of un- 
bounded non-helical flows with increasing rotation rate 
ft. As reported in p revious numerical simulations (2~H 
and experiments jl2l . [l3j , as il increases the decay slows 
down until reaching a saturated decay for Ro « 0.1. Wc 
will focus in the following in simulations with Rossby 
number small enough to observe this saturated decay, 
although not so small that the rotation quenches all non- 
linear interactions giving only exponential decay. A de- 
tailed study of the transition between the non-rotating 
and rotating cases can be found in [l2| . 
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FIG. 6: (a) Evolution of the isotropic energy spectrum E(k) 
for L96-6 (non-helical, Q = 10, initial ~ k 4 spectrum peaking 
at k = 14) from t = 5 to t = 100 with time increments At — 5. 
Inset: reduced perpendicular energy spectrum E(k±) for the 
same times, (b) Evolution of the isotropic energy spectrum 
for L96H-2 (helical, = 10, initial ~ k 4 spectrum peaking at 
k = 14) at the same times, with the reduced perpendicular 
energy spectrum in the inset. 



Figure [5] shows the energy, enstrophy, and helicity de- 
cay in simulations of rotating flows with and without 
helicity, in the unbounded and bounded cases (the latter 
in the insets). The energy decay in the unbounded non- 
helical runs (thin lines in Fig. [3]) is slightly steeper than 
what Eq. ([25]) predicts (E ~ i~ 5 / 7 ). A better agreement 
is observed for the enstrophy, which is closer to a ~ £~ 12 / 7 
law. As will be shown next, the agreement between the 
phcnomcnological arguments for the energy and the sim- 
ulations is improved if the decay of 2D and 3D modes is 
considered separately. 

Alternatively, the unbounded helical runs (thick lines 
in Fig. show for the energy a ~ i" 1 / 3 decay or steeper 
(although shallower than ~ t~ 5 / 7 ). Runs with ABC ini- 
tial conditions tend to develop a clearer power law de- 
cay and to be closer to a ~ i^ 1 / 3 decay than runs with 
random helical modes. Again, these differences can be 
explained considering the decay of 2D and 3D modes, as 
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FIG. 7: Axisymmetric energy spectrum e{k\\,k±)/ sm.0 for 
different times for run L96-6 (non-helical, Q — 10, initial en- 
ergy spectrum ~ k 4 peaking at k — 14). 



well as the effect of anisotropy in the initial conditions 
which is specially relevant for this particular case. The 
enstrophy and helicity show a decay close to ~ i~ 12 / 7 . 
Note that in the presence of rotation, helicity not only 
slows down the occurrence of the peak of enstrophy as 
already reported in |27j . but it also changes the energy 
decay after this peak. The enstrophy decay is not affected 
by the presence of helicity. 

Overall, the case of constrained runs shows a similar 
scenario, with a significant slow down of the decay rates 
in the presence of rotation, and with an extra slow down 
of energy decay in the presence of helicity (sec the in- 
sets in Fig. [3]). Rotating non- helical flows are close to 
E{t) ~ t-\ Q(i) ~ t~ 2 , and H(t) ~ t~ 2 , while heli- 
cal flows in this case display a shallower decay in the 
energy consistent with E ~ i -1 / 3 as predicted by the 
phenomenological arguments that take into account the 
effect of helicity in the energy spectrum of rotating t urbu- 
lence. As in the unbounded case, the presence of helicity 
does not affect the decay rate of enstrophy. 



2. Anisotropic global quantities 

Although the impact of rotation and helicity in the en- 
ergy decay is clear, the predictions given by the isotropic 
phenomenological arguments in Sec. lIIIBl do not coincide 
in all cases with the results from the simulations. This 
can be ascribed to the fact that these arguments assume 
an isotropic energy scaling, while rotation breaks down 
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FIG. 8: Axisymmetric energy spectrum e(fc||, k±)/sinO for 
different times for run L96H-3 (helical with ABC initial con- 
ditions, Q — 10, initial energy spectrum ~ k 4 peaking at 
k = 14). 



isotropy making spectral energy distribution become ax- 
isymmetric with the preferred direction along the axis 
of rotation. Two-dimcnsionalization of the flow has al- 
ready been reported during the decay 24|, 27 . as well 
as weak coupling of 2D and 3D modes [g, [l5( for very 
small Rossby numbers. Based on this, we discriminate 
between the energy contained in 3D modes with k z ^ 
[E^d], and the energy in slow 2D modes with k z = 
(E2d)- At this point it is important to point out that as 



the energy decays in the simulations, the Rossby number 
also monotonically decreases with time. As a rule, the 
Rossby numbers decrease by one order of magnitude in 
the first turnover times (from the values listed in Table HI 
which correspond to the time of maximum dissipation), 
and another order of magnitude at t ~ 100. 

In Fig. 2] we show E^d and E2D as a function of time for 
several runs. In each and every case we can clearly iden- 
tify distinct behaviors for the 2D and 3D energies, obey- 
ing different power-law decays. On the one hand, £3.0 
always shows a decay close to some power law expected 
for a (bounded or unbounded) non-rotating case, with 
the unbounded non- helical runs having E^d ~ £~ 10 / 7 
in agreement with Eq. (f20|) as illustrated in Fig. EJa), 
and with most helical runs (bounded and unbounded) 
with E-sd <~ t~ 2 in agreement with Eq. (fP9)) (which cor- 
responds to the bounded non-rotating decay) as illus- 
trated in Fig. @Jb-d). On the other hand, E^d follows 
power laws close to the ones predicted by Eqs. (pTLj) - 
(|37[) . The unbounded non- helical runs are compatible 



with E2D ~ t~ 2 / 3 , and the helical runs show ~ i -1 / 2 or 
~ £ -1 / 3 (note in the helical case the power laws predicted 
are for the case of maximum helicity, and for interme- 
diate helicity the power laws are bounded between the 
non- helical and the maximally helical values). 

The results in Fig. U indicate clearer power law de- 
cay (and better agreement with phenomcnological argu- 
ments) is obtained for the separate energy in 2D and 3D 
modes, than when the total energy is considered (com- 
pare, e.g., the extent of the power laws in this figure with 
the ones in Fig. [3]). This is clearer in the non- helical case, 
where all unbounded runs show a decay consistent with 
Eq. (f3Ul) for the 2D modes (all non-helical runs with ran- 
dom initial conditions have a <~ fc^ initial energy spec- 
trum, per virtue of the isotropic initial ~ fc 4 spectrum), 
and with Eq. (|20|) for the 3D modes. The separate evo- 
lution of E3D and E2D seems to be independent of the 
initial ratio of energy in 3D and 2D modes, at least for 
the range of parameters explored in this work. 

To further show the agreement with the phenomcno- 
logical arguments, the evolution of /, I2D and K must be 
considered, to verify whether these quantities behave in 
agreement with the assumptions used to derive the de- 
cay laws in Sec. lIIII Figure [5] shows the time evolution of 
I and I2D for two simulations with rotation, normalized 
by their initial values at t = 0. In both cases I grows 
monotonically in time by a factor of rj 50 during the de- 
cay. This fast increase of / is observed in all rotating 
runs. Meanwhile, I2D settles and remains approximately 
constant, fluctuating around its initial value. 

The helical runs show more disparity in the time evo- 
lution of the 2D and 3D energies. Bounded runs show 
E 2 d ~ i _1/3 and E 3D ~ i~ 2 , which agree with the 
previous scenario where 3D modes decay as in the non- 
rotating case, and slow 2D modes decay according to the 
anisotropic prediction with rotation (in this case, corre- 
sponding to Eqs. (|24p and (fTTj|) . respectively). But for 
initial conditions that peak at fc ps 14, in some cases they 
show decays of £"313 and E2D that are consistent with 
predictions for bounded flows (Fig. @] (b)), while in oth- 
ers they show decays as in the unbounded case (Fig. @] 
(c)). It may be the case that in the presence of helicity 
more separation of scales is needed between the initial 
energy containing scale and the size of the box in order 
to study unbounded flows (indeed, run L192H-2, which 
has an initial energy spectrum peaking at fc = 30, shows 
a decay compatible with E^jj ~ t~ 10 / 7 ). But it is also 
observed that these decay laws also depend on whether 
ABC or random helical initial conditions are used. In the 
ABC flow, two-thirds of the initially excited modes are in 
the fc|| = plane (see Eq. (|5J|), while random initial con- 
ditions excite modes in Fourier space distributed more 
isotropically, resulting in a smaller percentage of energy 
in the fc|| = modes when compared with the energy in 
fcii 7^ 0. This dependence in the initial ratio of energy 
in 2D and 3D modes may indicate a stronger coupling 
between 2D and 3D modes in the presence of helicity (in 
Sec. [V] we will explicitly show how changes in the initial 
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FIG. 9: Energy decay for different initial anisotropies. Thick 
lines correspond to the energy in 2D modes (E2d), while 
thin lines correspond to energy in 3D modes (E^d)- Simu- 
lations shown are L192HA-1 (solid), L192HA-2 (dotted), and 
L192HA-3 (dashed), with increasing anisotropy. 



anisotropy affect these results). 

On dimensional grounds, the impact of helicity in the 
coupling can be explained as follows. If decoupling takes 
place in the limit of fast rotation, it should hold until a 
time t* ~ Ro~ 2 , after which higher order terms in the 
multiple time scale expansion make non-resonant inter- 
actions relevant [HI, • In non- helical unbounded tur- 
bulence, E ~ i~ 5 / 7 and L ~ t 1 / 7 (under approximate 
conservation of I). The Rossby number then decays as 



Ro = 



E 1 ' 2 

¥J 2 m 



-1/2 



(39) 



and t* grows as t. As a result, if decoupling takes place in 
the freely decaying non-helical case, it can be sustained 
for long times. The same result (Ro ~ i^ 1 / 2 ) is obtained 
if the argument is refined to consider the 2D invariants 
K or I 2 d using Eqs. (|3"Tj) or (|33|) . or in the bounded case 
using Eq. (|2"2")l . However, in the helical case (e.g., using 
Eq. ([24)) ) a much slower decay of the Rossby number 
obtains 

Ro ~ r 1 / 6 , (40) 
and thus t* grows only as t 1 / 3 . 

C. Spectral evolution and anisotropy 

The isotropic and reduced perpendicular energy spec- 
tra are shown in Fig.|6]for LES of rotating flows (f2 = 10) 
with and without helicity. Energy at large scales grows 
in all cases, indicating an inverse energy transfer (as also 
evidenced by a negative flux of energy in that range). 
Also, the energy spectrum in the helical case (e.g., at 
the time of maximum enstrophy; not shown) is slightly 
steeper than in the non-helical case (see, e.g., (55l l5o|). 
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FIG. 10: Evolution of the velocity-derivative skewness for 
non-rotating runs (a) D512-2 and L96-1 (non-helical), and 
(b) D512H-2 and L96H-1 (helical). DNS have filled symbols 
while LES have empty symbols, with squares for S x , trian- 
gles for S y , and diamonds for S z . The three components of 
S oscillate around w —0.5 independently of helicity content. 
The insets show the three components of the kurtosis for the 
same runs using the same labels. 



To further investigate the energy spectral distribution 
among different directions, we show in Figs.[7]and[S]plots 
of the axisymmetric energy spectrum e(k\\,k±) for runs 
L96-6 and L96H-3. Note that to obtain contour levels 
corresponding to circles in the isotropic case, here and 
in the following, the axisymmetric energy spectrum is 
divided by sin 8, where 8 = arctan(fc|i/fcj_). 

In the case without rotation the spectrum has an 
isotropic distribution of energy evidenced by circular con- 
tour levels, which maintain their shape as the flow decays 
(not shown). Alternatively, when rotation is present, the 
distribution of energy becomes anisotropic with more en- 
ergy near the kn = axis at late times (Fig. [7]). This pref- 
erential transfer towards slow 2D modes is well known, 
see e.g., However, for helical rotating flows there 

is even more energy near the fcii = axis (Fig. [5]), and 
energy is also more concentrated at large scales (small k± 
wave numbers), in agreement with our previous observa- 
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tions of a faster increase of integral scales in the presence 
of helicity. 



V. EFFECT OF INITIAL ANISOTROPY 

As mentioned before, some of the differences observed 
in the evolution of E 2 d and E^u in helical runs are asso- 
ciated with differences in the initial conditions. In par- 
ticular, runs with ABC and with random helical initial 
conditions differ in the fact that the ABC flow initially 
has 2/3 of the excited modes in Fourier space in the slow 
2D manifold, while in the random case energy is more 
isotropically distributed. 

To further investigate this effect, we consider a set of 
helical runs with random initial conditions but with in- 
creasing initial anisotropy (runs L192HA-1, L192HA-2, 
and L192HA-3). Anisotropy is introduced by weighting 
the amplitude of all modes with fen = with a parame- 
ter a (a = 1 corresponds to the isotropic initial condi- 
tions considered before, and a > 1 corresponds to larger 
amplitude of the 2D modes relative to the 3D modes). 
The runs have a = 1 (L192HA-1), 5 (L192HA-2), and 
10 (L192HA-3), resulting in initial ratios of the energy 
in 2D to 3D modes E 2D /E 3D ps 0.024, 0.626, and 2.408, 
respectively. The runs also have initial energy and helic- 
ity spectra peaking at k = 25, thus allowing us to study 
unbounded cases with larger scale separation. 

Figure [5] shows that E^o decays approximately as 
ps —10/7 regardless of the anisotropy of the initial con- 
ditions, while E2D changes its decay becoming shallower 
as anisotropy grows. On the one hand, the isotropic case 
(a = 1) is closer to a E 2 d ~ t~ x l z law, a result consis- 
tent with the 2D decay shown in Fig. HJd) . On the other 
hand, the decay in the most anisotropic case (a = 10) 
is closer to ~ i -1 / 6 , which is consistent with the decay 
for helical flows in the case when K is approximately 
constant; see Eq. (f53|) . Indeed, it was verified that K 
remains approximately constant during the decay of this 
run (not shown). 



VI. SKEWNESS AND KURTOSIS 

In this section we consider the time evolution of skew- 
ness and kurtosis of velocity derivatives in runs with and 
without rotation, and with and without helicity. The 
skweness St and kurtosis Ki are defined as 
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(41) 




(42) 



1.0 - 



0.5 



3.6 



K gKB h h 



[a) _ 



20 40 60 



100 



A A A A A 



20 



40 60 80 100 

time 



2 r 



5 : 

4 
2 




b z 



AAA. 



40 60 
time 

n □ □ 



100 

A - 



A 



AA n 



20 



40 60 
time 



80 



100 



where i denotes the Cartesian coordinates x, y, or z. 



FIG. 11: Evolution of the velocity-derivative skewness for 
runs (a) L96-6 and (b) L96H-3. Symbols are squares for S x , 
triangles for S y , and diamonds for S z . The inset shows the 
evolution of velocity-derivative kurtosis for the same runs. 



Figure [TOl shows S and K for non-rotating runs D512-2 
and L96-1 (non-helical), and D512H-2 and L96H-1 (heli- 
cal). Only a few times are shown for the DNS runs, to 
compare with the LES. Overall, the DNS and LES show 
good agreement. The three components of the skewness 
fluctuate around ~ —0.5, a value observed in experiments 
[57] | and simulations [58(. Also, the kurtosis evolves to- 
wards a value near 3.5. Helicity does not affect the values 
of S nor K in the absence of rotation. 

When rotation is present skewness is substantially 
reduced, with all components of S fluctuating around 
S ps 0. This is shown in Fig. [11] for runs L96-6 and 
L96H-3 (DNS show the same behavior and are not shown 
for clarity). Such a decrease of skewness with decreasing 
Rossby number has already been reported in simulations 
Q. Anisotropy is also evident, manifested in a distinct 
behavior of S x , S y , and S z . While fluctuations of S z are 
small, S x and S y show large and sudden departures from 
zero with S x ps —S y at all times. This anti-correlation 
between the x and y components can be qualitatively 
understood from the two-dimcnsionalization of the flow. 
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FIG. 12: (Color online) Visualization of uj z at late times for 
run L96H-3. From top to bottom and from left to right, the 
images correspond to t = 42.5, 47.5, 55, 70, t = 87.5, and 100. 
Note that the four anti-cyclonic vortices at t — 42.5 merge in 
pairs and become two larger columnar vortices at t = 47.5. 
Eventually they merge again becoming one column. 



For a 2D flow, the incompressibility condition becomes 



du x 
dx 



du 



y_ 

dy 1 



(43) 



which leads to S x ~ —S y . 

The kurtosis in the runs with rotation has more fluc- 
tuations, but seems to evolves towards a value near 3. 
This is clearer for K z , while K x and K y also show signs 
of two-dimensionalization with K x « K y at all times. 

Visualization of the flow vorticity indicates that max- 
ima and minima of S x and S y correspond to times when 
two column-like structures in the flow merge. As an ex- 
ample, Fig. U2l shows the evolution of the z component of 
the vorticity in run L96H-3. When two columns with the 
same sign of vorticity merge, intense gradients are cre- 
ated in u x and u y , giving rise to an increase or decrease 
in the values of S x and S y . Columns of positive vortic- 
ity (cyclonic) tend to merge, while vortices of negative 
vorticity (anti-cyclonic) seem to be unstable. 



VII. CONCLUSIONS 

In this work we presented a study of the self-similar de- 
cay laws that arise in turbulent rotating flows depending 
on: (1) the characteristic scale of the initial conditions 
(compared with the size of the box), (2), the presence 
or absence of helicity in the flow, (3) the values of the 
Rossby and Reynolds numbers, and (4) the amount of 
anisotropy in the initial conditions. Phenomenological 
decay laws were obtained for each case considered, and 
the decay laws were contrasted with numerical results 
from DNS and LES using different flows as initial condi- 
tions. 

A large number of power laws were identified. It is 
well known that rotation decreases the energy decay rate 
0, US EHH, HI HH , and our simulations are in agree- 
ment with this result. However, our simulations fur- 
ther show that in the presence of rotation helicity can 
further decrease this decay. This is different from the 
non-rotating case, where helicity does not affect the self- 
similar decay of energy. This result, together with previ- 
ous studies in the case of forced rotating flows [ll], [Ffl [56| 
further confirm that helicity plays a more important role 
in rotating turbulence than what it docs in the isotropic 
and homogeneous case. 

In the presence of rotation, the decay of enstrophy is 
well described by phenomenological arguments based on 
isotropic scaling. This can be expected as enstrophy (as 
well as helicity) is a small-scale quantity, more isotropic 
than the energy. 

The energy in rotating non-helical flows follows either 
a decay close to a ~ t~ 1 law (when the integral scale 
of the flow is close to the size of the box), or a decay 
slightly steeper than ~ t™ 5 / 7 (when the integral scale 
is smaller than the size of the box, and the large scale 
energy spectrum is ~ fc 4 ). Better agreement with power- 
law decay is obtained when the evolution of 2D modes 
and 3D modes is considered separately. In that case, the 
energy in 2D modes decays close to E 2 d ~ t~ 2 / 3 , and the 
3D modes decay as in the non- rotating case, i.e., close to 
E 3D ~ t~ 10 /7. 

These power-law decays can be obtained from phe- 
nomenological arguments that consider the energy in 2D 
and 3D modes separately, that assume approximately 
constant axisymmetric integrals instead of the isotropic 
Loitsyanski's integral, and that take into account the 
slow down in the energy transfer associated with rota- 
tion. Note we are not claiming there is decoupling be- 
tween 2D and 3D modes in rotating flows, a topic which 
is beyond the scope of this work. What we show instead 
is that the energy in 2D and 3D modes in the simulations 
decay with different rates, both following power laws, and 
that considering this in the phenomenological description 
gives better agreement with the numerical results. 

The decay of energy in the presence of rotation and 
helicity shows further variety. When the integral scale of 
the flow is close to the size of the box, the energy decay 
is close to E ~ tr 1 ^ 3 . This decay can be obtained from 
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phenomenological arguments taking into account the role 
played by the helicity cascade in further slowing down the 
energy transfer. In simulations with integral scale smaller 
than the size of the box, the decay is between E ~ t -5 / 7 
and ~ t -1 / 3 . As in the non-helical case, clearer power 
laws arise if the decay of E^d and E^d is considered. 
In that case, E^d shows decays between ~ t^ 1 / 2 and 
~ t -1 / 6 , and E3D shows decays close to either E ~ t~ 2 
or E ~ t~ w / 7 . 

The results with helicity seem to be more dependent 
on scale separation (i.e., on the separation between the 
initial integral scale of the flow, and the size of the box), 
and on initial anisotropy. It is worth mentioning that 
the importance of the initial conditions in the decay of 
rotating turbulence has been recently pointed out also in 
experiments [59l | . 

Finally, we presented a study of the time evolution of 



the skewness and kurtosis of velocity derivatives. Two- 
dimcnsionalization of rotating flows leads to an anti- 
correlation of the x and y components of the skewness, 
which fluctuate around zero. Large departures of these 
quantities from this value are associated with merging 
events of columns in the flow. 
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